#include "HY1C_out.h"
#include "l12_proto.h"

#include <math.h>

/* ========================================================================== */
/* module atmocor1() - computes pre-aerosol atmospheric components            */
/*                                                                            */
/* Written By: B. A. Franz, SAIC, SeaWiFS Project, February 1998              */
/* Conversion to C, May 2006                                                  */
/* ========================================================================== */


/* ------------------------------------------------------------------- */
/* correction factor to remove oxygen absorption from La(765)          */
/* ------------------------------------------------------------------- */
float oxygen_aer(float airmass)
{
    /* base case: m80_t50_strato  visibility (550nm ,0-2km):25km */
    static float a[] = {-1.0796,9.0481e-2,-6.8452e-3};
    return(1.0+pow(10.0,a[0]+airmass*a[1]+airmass*airmass*a[2]));
}


/* ------------------------------------------------------------------- */
/* correction factor to replace oxygen absorption to Lr(765)           */
/* ------------------------------------------------------------------- */
float oxygen_ray(float airmass)
{
    /* base case is the 1976 Standard atmosphere without aerosols */
    static float a[] = {-1.3491, 0.1155, -7.0218e-3};
    return(1.0/(1.0+pow(10.0,a[0]+airmass*a[1]+airmass*airmass*a[2])));
}	


#ifndef PI 3.14159265358979323846
#define PI 3.14159265358979323846
#endif

static float HY1C_acoss(float x1){
    float a;
    if (x1>=1.0){
        a=0.0;
    }
    else{
        if (x1<=-1.0){
            a=3.141592654;
        }
        else{
            a=acos(x1);
        }
    }
    return a;
}

static void HY1C_reflec_both(float x1, float *out1, float *out2){
    float x2;
    float x3;
    float x4;
    float  ref=4.0/3.0;
    if (x1<0.00001){
        x3=0.0204078;
        x4=0.0;
    }
    else{
        float x2=asin(sin(x1)/ref);
        float perp=pow((sin(x1-x2)/sin(x1+x2)),2);
        float par=pow((tan(x1-x2)/tan(x1+x2)),2);
        x3=perp+par;
        x4=-perp+par;
        x4=x4/2.0;
    }
    //one_array result(2,{x3,x4});
    *out1 = x3;
    *out2 = x4;
    //return result;
}

static float HY1C_asinn(float x1){
    float a;
    if (x1>=1.0){
        a=3.141592654/2.0;
    }
    else{
        if (x1<=-1.0){
            a=-3.141592654/2.0;
        }
        else{
            a=asin(x1);
        }
    }
    return a;
}

static void HY1C_get_glint(float x1, float x2, float x3, float x4, float x5, float *out1, float *out2, float *out3){
    //one_array resultres(2);
    float y4=fmax(x4,0.001);
    float y1=x1*PI/180.0;
    if (y1==0.0)
        y1=1.0e-7;
    float y2=x2*PI/180.0;
    if (y2==0.0)
        y2=1.0e-7;
    float y3=x3*PI/180.0;
    float omega1=HY1C_acoss(cos(y1)*cos(y2)-sin(y1)*sin(y2)*cos(y3))/2.0;
    if (omega1==0.0)
        omega1=1.0e-7;
    float beta1=HY1C_acoss((cos(y1)+cos(y2))/(2.0*cos(omega1)));
    if (beta1==0.0)
        beta1=1.0e-7;
    float alpha1=HY1C_acoss((cos(beta1)*cos(y2)-cos(omega1))/(sin(beta1)*sin(y2)));
    if (sin(y3)<0.0)
        alpha1=-1.0*alpha1;
    float sigc=0.04964*sqrt(y4);
    float sigu=0.04964*sqrt(y4);
    float chi=x5;
    float alphap=alpha1+chi;
    float swig=sin(alpha1)*tan(beta1)/sigc;
    float eta=cos(alpha1)*tan(beta1)/sigu;
    float expon=-(pow(swig,2)+pow(eta,2))/2.0;
    if (expon<-30.0)
        expon=-30.0;
    if (expon>30.0)
        expon=30.0;
    float prob=exp(expon)/(2.0*PI*sigu*sigc);
    //resultres=reflec_both(omega1);
    float rho_plus = 0;// resultres.data[0];
    float rho_minus = 0;// resultres.data[1];
    HY1C_reflec_both(omega1, &rho_plus, &rho_minus);
    float x6=rho_plus*prob/(4.0*cos(y1)*(pow(cos(beta1),4)));
    float rot_ang;
    if (omega1>0.0001){
        float cr=(cos(y2)-cos(2.0*omega1)*cos(y1))/(sin(2.0*omega1)*sin(y2));
        float sr=sin(y2)*sin(PI-y3)/sin(2.0*omega1);
        rot_ang=(cr/fabs(cr))*HY1C_asinn(sr);
    }
    else{
        rot_ang=PI/2.0;
    }
    float c2r=cos(2.0*rot_ang);
    float s2r=sin(2.0*rot_ang);
    float x7=c2r*rho_minus/rho_plus;
    float x8=-s2r*rho_minus/rho_plus;

    // one_array result(3,{x6,x7,x8});

    *out1 = x6;
    *out2 = x7;
    *out3 = x8;

    //return result;
}

/* ------------------------------------------------------------------- */
/* main function, loads various quanitities into l1 record for 1 pixel */
/* ------------------------------------------------------------------- */

void atmocor1(l1str *l1rec, int32_t ip)
{
    static float  p0 = STDPR;
    static float  pi = PI;
    int32_t  sensorID = l1rec->sensorID;
    // 20180822
    int32_t  nwave    = 8;//l1rec->nbands;

    float solz = l1rec->solz[ip];
    float senz = l1rec->senz[ip];
    float raz  = l1rec->delphi[ip];
    float mu0  = l1rec->csolz[ip];
    float mu   = l1rec->csenz[ip];

    float ws   = l1rec->ws[ip];
    float pr   = l1rec->pr[ip];
    float oz   = l1rec->oz[ip];
    float wv   = l1rec->wv[ip];
    int nbands=8;
    float no2_tropo = l1rec->no2_tropo[ip];
    float no2_strat = l1rec->no2_strat[ip];
    float no2_frac = l1rec->no2_frac[ip];

    int32_t  *wave  = l1rec->iwave;
    float *Fo    = l1rec->Fo;
    float *Tau_r = l1rec->Tau_r;

    instr *input = l1rec->input;

    float zero = 0.0;
    int32_t  ib765 = -1;
    float airmass;
    float a_o2;
    float glint_coef_q;
    float glint_coef_u;
    int32_t  ib, ipb;

    airmass = 1.0/mu0 + 1.0/mu;

    /* Initialize output values */

    for (ib=0; ib<nwave; ib++) {

        ipb = ip*nbands + ib;

        l1rec->Lr  [ipb] = 0.0;
        l1rec->L_q [ipb] = 0.0;
        l1rec->L_u [ipb] = 0.0;
        l1rec->tLf [ipb] = 0.0;
        l1rec->tg_sol[ipb] = 1.0;
        l1rec->tg_sen[ipb] = 1.0;
        l1rec->t_h2o   [ipb] = 1.0;
        l1rec->t_o2    [ipb] = 1.0;
        l1rec->t_sol   [ipb] = exp(-0.5*pr/p0*Tau_r[ib]/mu0);
//        float check1=l1rec->t_sol   [ipb];
        l1rec->t_sen   [ipb] = exp(-0.5*pr/p0*Tau_r[ib]/mu );
//        float check2=l1rec->t_sen   [ipb];
        if (input->oxaband_opt > 0 && wave[ib] == 765) {
	    ib765 = ib;
            l1rec->t_o2[ipb] = 1.0/oxygen_aer(airmass);
	}

        if (sensorID == SEAWIFS) {
	    /* this is for rhos only, but effects seawifs land products and cloud   */
            /* will modify get_rhos() to use gaseous transmittance in future update */
	    l1rec->t_h2o[ipb] = water_vapor(ib,wv,airmass);
        }

    }

    /* gaseous transmittance */

    ipb = ip*nbands;
    gaseous_transmittance(input->gas_opt,sensorID,input->evalmask,nwave,mu0,mu,oz,wv,no2_tropo,no2_strat, no2_frac,
                          &l1rec->tg_sol[ipb],&l1rec->tg_sen[ipb]);

    /* white-cap radiances at TOA */

    ipb = ip*nbands;
    whitecaps(sensorID,input->evalmask,nwave,ws,input->wsmax,&l1rec->rhof[ipb]);
    for (ib=0; ib<nwave; ib++) {
        ipb = ip*nbands+ib;
        l1rec->tLf[ipb] = l1rec->rhof[ipb]*l1rec->t_sen[ipb]*l1rec->t_sol[ipb]*Fo[ib]*mu0/pi;

    }

    /* Rayleigh scattering */

    ipb = ip*nbands;
//    sd_id = SDstart("/home/oceanoptics/seadas64/run/data/modisa/rayleigh/rayleigh_modisa_412_iqu.hdf", DFACC_RDONLY);
//HY1C_out("%d",sd_id);
    rayleigh(sensorID,input->evalmask,nwave,input->pol_opt,solz,senz,raz,Tau_r,Fo,pr,ws,
             &l1rec->Lr[ipb],&l1rec->L_q[ipb],&l1rec->L_u[ipb]);             

    if (input->oxaband_opt > 0 && ib765 > -1) {
        a_o2 = oxygen_ray(airmass);
        l1rec->Lr [ipb+ib765] *= a_o2;
        l1rec->L_q[ipb+ib765] *= a_o2;
        l1rec->L_u[ipb+ib765] *= a_o2;
    }


    /* glint coefficients and approximate glint radiances */
    /* also add glint to polarization components          */

    // getglint_iqu_(&senz,&solz,&raz,&ws,&zero, &l1rec->glint_coef[ip],&glint_coef_q,&glint_coef_u);
    HY1C_get_glint(senz,solz,raz,ws,zero, &l1rec->glint_coef[ip],&glint_coef_q,&glint_coef_u);

    for (ib=0; ib<nwave; ib++) {
        ipb = ip*nbands+ib;
//        l1rec->glint_coef[ip]=0.1;
        l1rec->TLg[ipb]  = l1rec->glint_coef[ip] * exp(-(Tau_r[ib]+0.1)*airmass)*Fo[ib];
        l1rec->L_q[ipb] += (glint_coef_q * l1rec->TLg[ipb]);
        l1rec->L_u[ipb] += (glint_coef_u * l1rec->TLg[ipb]);
    }

}


